function Sumrate = Sumrate_Detequi(K,L,N,T2,Hbar2,tau,rho,alpha,H1,F)
for ueIdx = 1:K
    Hbar(ueIdx,:) =  Hbar2(ueIdx,1:L) * F * H1;
    T(:,:,ueIdx) =  H1' * F' * T2(:,:,ueIdx) * F * H1; 
    %SQR_T(:,:,ueIdx) =  sqrtm(T(:,:,ueIdx));
end
[e,Psi,tPsi] = Itrative_DeterEq(K,L,N,T,Hbar,alpha);

u0 = trace(Psi)/L;
v0 = vbar(eye(N),eye(N),K,L,N,T,Hbar,e,Psi,tPsi,alpha);  

D_noise = u0 - alpha*v0;

for k = 1:K
    u1(k) = trace(Hbar(k,:)'*Hbar(k,:)*Psi);
    u2(k) = trace(T(:,:,k)*Psi)/L;       
    v1(k) = vbar(Hbar(k,:)'*Hbar(k,:)*L,eye(N),K,L,N,T,Hbar,e,Psi,tPsi,alpha);
    v2(k) = vbar(T(:,:,k),eye(N),K,L,N,T,Hbar,e,Psi,tPsi,alpha);
    
    a(k) = (u1(k)+sqrt(1-tau(k)^2)*u2(k))/(1+u2(k)) - sqrt(1-tau(k)^2)*u1(k)*u2(k)/(1+u2(k))^2;
    b(k) = ((1+u2(k))^2*(u1(k)-alpha*v1(k))-u1(k)^2*(1+u2(k)-alpha*v2(k)))*(a(k)-1)^2/(1+u2(k)-u1(k))^2 ...
        + (u2(k)-alpha*v2(k))*(a(k)^2-2*sqrt(1-tau(k)^2)*a(k)+1);
       
    D_signal(k) = rho*a(k)^2;
    D_interferece(k) = rho*b(k);
    
%     for kk = 1:K
%         v3(k,kk) = vbar(Hbar(kk,1:Nt)'*Hbar(kk,1:Nt)*Nt,Hbar(k,1:Nt)'*Hbar(k,1:Nt)*Nt,K,Nt,T,Hbar,e,Psi,tPsi,alpha);
%         v4(k,kk) = vbar(T(1:Nt,1:Nt,kk),Hbar(k,1:Nt)'*Hbar(k,1:Nt)*Nt,K,Nt,T,Hbar,e,Psi,tPsi,alpha);
%         v5(k,kk) = vbar(Hbar(kk,1:Nt)'*Hbar(kk,1:Nt)*Nt,T(1:Nt,1:Nt,k),K,Nt,T,Hbar,e,Psi,tPsi,alpha);
%         v6(k,kk) = vbar(T(1:Nt,1:Nt,kk),T(1:Nt,1:Nt,k),K,Nt,T,Hbar,e,Psi,tPsi,alpha);
%     end     
 end
% 
% for k = 1:K
%     for kk = 1:K
%         tempu3(k,kk) = (a(kk)-1)^2*((1+u2(kk))^2*v3(k,kk)-u1(kk)^2*v4(k,kk))/(1+u2(kk)-u1(kk))^2 + (a(kk)^2-2*sqrt(1-tau(kk)^2)*a(kk)+1)*v4(k,kk);
%         tempu4(k,kk) = (a(kk)-1)^2*((1+u2(kk))^2*v5(k,kk)-u1(kk)^2*v6(k,kk))/(1+u2(kk)-u1(kk))^2 + (a(kk)^2-2*sqrt(1-tau(kk)^2)*a(kk)+1)*v6(k,kk);
%     end
% end
% for k = 1:K
%     u3(k) = (sum(tempu3(k,:)) - tempu3(k,k))/Nt;
%     u4(k) = (sum(tempu4(k,:)) - tempu4(k,k))/Nt;  
% end

for k = 1:K
    De_gamma1(k) = D_signal(k)/(D_interferece(k)+D_noise);
    temp_rate1(k) = log2(abs(1+De_gamma1(k)));
    
%     D_eavesd(k) = rho*((u3(k)+u4(k))*(1+u2(k))-2*u1(k)*u4(k))/(1+u2(k))^3;
%     
%     De_gamma2(k) = D_eavesd(k)/D_noise;
%     temp_rate2(k) = log2(abs(1+De_gamma2(k)));
%     
%     temp_Secreyrate(k) = temp_rate1(k) - temp_rate2(k);
%     if temp_Secreyrate(k)<=0
%         temp_Secreyrate(k) = 0;
%    end
end
Sumrate =  sum( temp_rate1);

%TempD = [abs(D_noise),abs(D_signal(3)),abs(D_interferece(3)),abs(D_eavesd(3))]    %user 3
end
%% ---------------Itrative ----------------
function [e,Psi,tPsi] = Itrative_DeterEq(K,L,N,T,Hbar,alpha)

%---------Initialize----------
e = ones(K,1);
tPsi = eye(K);
%---------Iterative-----------
t=1;
old_e = e+ones(K,1)*1000;
Psi = zeros(N,N);  %*
while (sum(sum(abs(old_e - e))) > 10^(-8))
   old_e = e; 
   t=t+1;
   
   temp_Phi = eye(N);
   for k =1:K
       temp_Phi = temp_Phi + tPsi(k,k)*T(:,:,k)/L;
   end
   Psi = inv(alpha*temp_Phi+Hbar'*diag(1./(e+1))*Hbar);
   tPsi = inv(alpha*eye(K)+alpha*diag(e)+Hbar*inv(temp_Phi)*Hbar');
   
   for k = 1:K 
       e(k) = trace(T(:,:,k)*Psi)/L;   
   end
end
end

%% -----------------v Function--(Lemma 1)-----------------------
function v = vbar(Q1,Q2,K,L,N,T,Hbar,e,Psi,tPsi,alpha)
temp_phi = eye(N); 
temp_matrix1 = zeros(K,K,K);
temp_matrix2 = zeros(K,K,K);

for k = 1:K
   temp_phi = temp_phi + tPsi(k,k)*T(:,:,k)/L;
end

for k = 1:K
   temp_matrix1(:,:,k) = Hbar*Psi*T(:,:,k)*Psi*Hbar';
   temp_matrix2(:,:,k) = tPsi*Hbar*inv(temp_phi)*T(:,:,k)*inv(temp_phi)*Hbar'*tPsi; 
end
temp_matrix3 = tPsi*Hbar*inv(temp_phi)*Q2*inv(temp_phi)*Hbar'*tPsi; 

for k = 1:K
    temp_epsilon1(k,1) = trace(T(:,:,k)*Psi*Q2*Psi)/L;
    temp_epsilon2(k,1) = -temp_matrix3(k,k)/alpha;
    for kk = 1:K
        temp_Gamma11(k,kk) = temp_matrix1(kk,kk,k)/(1+e(kk))^2/L;
        temp_Gamma12(k,kk) = alpha*trace(T(:,:,k)*Psi*T(:,:,kk)*Psi)/L^2;
        temp_Gamma21(k,kk) = alpha*tPsi(k,kk)*tPsi(kk,k);
        temp_Gamma22(k,kk) = temp_matrix2(k,k,kk)/L;
    end
end

dotC = inv([eye(K)-temp_Gamma11,temp_Gamma12;temp_Gamma21,eye(K)-temp_Gamma22])*[temp_epsilon1;temp_epsilon2];

temp_num1 = 0;
temp_num2 = 0;

temp_matrix4 = Hbar*Psi*Q1*Psi*Hbar';

for k = 1:K
    temp_num1 = temp_num1 + temp_matrix4(k,k)*dotC(k)/(1+e(k))^2; 
    temp_num2 = temp_num2 + trace(Q1*Psi*T(:,:,k)*Psi)*dotC(K+k); 
end
v = trace(Q1*Psi*Q2*Psi)/L + temp_num1/L - temp_num2*alpha/L^2;
end

